//
// This file is part of cpp-lib, a library of C++ functions and classes.
// Cpp-lib is copyright (C) 2004 and onwards, The CPL Interest Group.
// 
// This library is free software; you can redistribute it and/or modify it
// under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation; either version 2.1 of the License, or (at
// your option) any later version.
// 
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
// General Public License for more details.
// 
// You should have received a copy of the GNU Lesser General Public License
// along with this library (see the file COPYING); if not, write to the
// Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
// 02111-1307  USA
//
// $Id$
//


#ifndef CPP_LIB_RANDOM_H
#define CPP_LIB_RANDOM_H

#include <cstdlib>
#include <cmath>
#include <cassert>

// Some rudimentary Pseudo Random Number Generation.


namespace cpl {
namespace math {

struct system_rng {

  double operator()() // NOT const!!!
  { return std::rand() / ( RAND_MAX + 1. ) ; }

} ;


/// Return X ~ lambda * exp( -lambda x )

template< typename rng >
double exponential_distribution( rng& r , double const& lambda ) {

  assert( lambda > 0 ) ;
  double const u = 1 - r() ;
  assert( u >  0 ) ;
  assert( u <= 1 ) ;

  return -std::log( u ) / lambda ;

}



//
// Return
//
//   X = sigma * ( X_1 + ... + X_n - n / 2 ) / sqrt( n / 12 ),
//
// where X_i ~ U[ 0 , 1 ] and independent.
//
// Note that for U[0,1], sigma^2 = 1/12, hence for X_1 + ... + X_n it's n/12.
//
// By the Central Limit Theorem, X's distribution will approach
// N( 0 , sigma ) as \a n -> infinity.
//

template< typename rng >
double n_times_distribution( rng& r , double const& sigma , int const n ) {

  assert( n >= 1 ) ;

  double sum = 0 ;
  for( int i = 0 ; i < n ; ++i ) { sum += r() ; }

  return sigma * ( sum - n / 2. ) / std::sqrt( n / 12. ) ;

}

} // end namespace math

} // end namespace cpl


#endif // CPP_LIB_RANDOM_H
